/*
    Copyright 2006-2008 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

// $Id$


#include "interpolatort.h"
#include "lsqfitter2.h"
#include "mcrx-types.h"

using namespace blitz;
using namespace mcrx;

typedef interpolator< T_float , T_float , 1> T_interpolator;


class points: public lsqfitter2 {
public:
  T_interpolator big,i;
  points (): lsqfitter2 (1e-5) {
    setFailretries (20);
    setAccuracy (1e-5,1+1e-7 );
    setSuccessincr (.5);
  }
  double function (double x, const matr& xpoints);
};

double points::function (double x, const matr& xpoints)
{
  vector<T_float> xvector;
  for (int j = 0; j < xpoints.getrows(); ++j) {
    const T_float current = xpoints [j];
    xvector.push_back(current );
  }

  i.setAxis(0, xvector );
  for (int j = 0; j < xvector.size(); ++j)
    i.setPoint(j, big.interpolate(xvector [j] ));
  
  return i.  interpolate (x);
}

int main (int argc, char**argv)
{
  // read file
  ifstream file (argv [1]);
  assert (file);
  vector<T_float> wavelengths ((istream_iterator<T_float> (file)),
			       istream_iterator<T_float> ());
  file.close ();
  file .clear();
  file.open(argv [2]);
  assert (file);
  vector<T_float> attenuation ((istream_iterator<T_float> (file)),
			       istream_iterator<T_float> ());
  file.close ();
  file .clear();
  file.open(argv [3]);
  vector<T_float> initial ((istream_iterator<T_float> (file)),
			       istream_iterator<T_float> ());
  file.close ();

  assert (wavelengths.size() ==attenuation.size());
  assert (wavelengths.size());
  
  for (int i = 0; i < wavelengths.size(); ++i) {
    wavelengths [i] = log10 (wavelengths [i]);
    attenuation [i] = log10 (attenuation [i]+ blitz::tiny (T_float ()));
  }
  for (int i = 0; i <initial.size(); ++i)
    initial [i] = log10 (initial[i]);
  
  points p;

  ofstream out ("fl");

  // the interpolator that describes the "true" function
  const T_float minimum = log10 (2.1e-8), maximum = log10 (4.98e-6);
  const int n_big = 1000;
  p.big.setAxis(0, wavelengths);
  for (int i = 0; i < attenuation.size(); ++i) {
    out << wavelengths [i] << "\t" << attenuation [i] << endl;
    p.big.setPoint(i, attenuation [i]);
  }
  arr xvalues(n_big), yvalues (n_big), sigma (n_big);
  int i=0;
  for (T_float x = minimum; i < n_big;
       x+= (maximum -minimum)/(n_big-1), ++i) {
    xvalues [i] = x ;
    yvalues [i] = p.big.interpolate(x);
    sigma [i] = .01;

    //cout << xvalues [i] << endl;
  }
  
  // generate initial points
  int n;
  if(initial.size())
    n=initial.size();
  else {
    sscanf(argv [3], "%i", & n);
  }
  arr xp (n);
  table<char>  fix (n);
  if(initial.size())
    for(int i=0;i<n;++i)
      xp[i]=initial[i];
  else {
    i=0;
    for (T_float x = minimum; i < n;
	 x+= (maximum -minimum)/(n-1), ++i) {
      xp [i] = x;
    }}
  fix=1;
  xp [0] = minimum -.001;
  xp [n- 1] = maximum+.001;
  fix [0] = 0;
  fix [n- 1] = 0;

  cout << n << endl;
  out << endl << endl;
  for (int i = 0; i < n; ++i) {
    out << xp [i] << " " << p.big.interpolate(xp [i]) << endl;
  }

  p.fit(xvalues, yvalues, sigma, xp, fix);
  out << endl << endl;
  for (int i = 0; i < n; ++i) {
    out << xp [i] << " " << p.big.interpolate(xp [i]) << endl;
  }
  
}
